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1 Abstract 
On 

Path integrals are usually formulated in discrete Euclidean time using the 

B 

i . Trotter formula. We propose a new method to study discrete quantum sys- 

terns, in which we work directly in the Euclidean time continuum. The method 

O ■ 

is of general interest and can be applied to studies of quantum spin systems, 
^ . lattice fermions, and in principle also lattice gauge theories. Here it is applied 

to the Heisenberg quantum antiferromagnet using a continuous-time version 
of a loop cluster algorithm. The computational advantage of this algorithm 
is exploited to confirm the predictions of chiral perturbation theory in the 
extreme low temperature regime, down to T = 0.01J. A fit of the low-energy 
parameters of chiral perturbation theory gives excellent agreement with pre- 
vious results and with experiments. 
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The path integral for a general quantum system is usually formulated in discrete Eu- 
clidean time. Farhi and Gutmann pointed out that this is not necessary, since the path 
integral is well-defined in continuous time [[[]] . For example, for a non-relativistic particle in 
coordinate space, each trajectory in the path integral resembles that of a particle undergoing 
Brownian motion. Due to the continuum of position eigenstates, defining the corresponding 
Wiener measure seems to require discrete time. However, one can always work in a count- 
able basis of the Hilbert space: one can put the particle in a large box and work in the 
discrete basis of momentum eigenstates, for example. The important point is that then the 
trajectories in the path integral do not resemble Brownian motion. Instead they consist of 
segments for which the system is in a basis state for a finite time, with sporadic discrete 
transitions between basis states. Consequently, for discrete quantum systems the path inte- 
gral can be formulated directly in continuous time. This insight can be applied, for example, 
to quantum spin systems, lattice fermions, and lattice gauge theories with a compact gauge 
group. 

In the continuous-time formulation a path is characterized by transition times and in- 
formation about which states are connected at the transitions. This picture allows the path 
integral to be sampled numerically without having to store information about the individual 
time slices whose number diverges in the time-continuum limit. In the discrete-time approx- 
imation one must always exercise care in extrapolating the results to the time-continuum 
limit. The approach advocated here cuts this Gordian knot by operating directly in this 
limit. This completely eliminates the most severe systematic error in these calculations. 

As an example we consider quantum spin systems. Conventional approaches to handling 
these systems rely on a discrete-time formalism. Consider, as a concrete example, the 
Hamiltonian for the Heisenberg quantum antiferromagnet 

h = jy^s x ■ s x+ ft, (i) 

where S x = \a x is a spin 1/2 operator associated with the site x of a d- dimensional hy- 
percubic lattice with spacing a. The interaction is between nearest neighbors; ft is the unit 
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vector in the //-direction. For an antiferromagnet the coupling constant J is positive. Since 
this Hamiltonian comprises non-commuting terms, the explicit evaluation of the associated 
partition function is problematical. Suzuki is credited with showing how the Trotter formula 
can be applied to segregate these non-commuting terms into separate time slices, with the 
approximation becoming exact in the continuum limit [TJ]. For example, the Hamiltonian of 
a ID quantum spin chain can be decomposed into two terms, H — H± + H2, each of which 
comprises only commuting terms 

Hi = J S x - S x+ i, H 2 = J ' $x+i- ( 2 ) 

x=2m x=2m+l 

Then for the partition function one writes 

Z = Trexp (-/?#) 

= lim Trlexpi-efiH^ exp(-e[3H 2 )} N (3) 

iv— »oo 

where (3 = 1/T is the inverse temperature and e/3 = f3/N is the lattice spacing in Eu- 
clidean time. Inserting complete sets of eigenstates |1) and | — 1) of between the factors 
exp(—e(3Hi) one converts the partition function into a (d + l)-dimensional path integral of 
Ising-like variables 

Z = H J2 «p(-S). (4) 

x,t s(x,t)=±l 

The action is a sum of plaquette couplings 

S = S[s(x,t),s(x + l),s(x,t+i),s(x + i,t + l)] 

x=2m,t=2p 

+ Y S[s{x,t),s(x + l),s(x,t + l),s{x + i,t + l)], 

x=2m+l,t=2p+l 

(5) 

and the plaquette Boltzmann factors are given by a 4 x 4 transfer matrix 

exp(-S , [si,S2,s 3 ,s 4 ]) = (sis 2 | exp(-e(3JS x -S x+fl )\s 3 S4). (6) 

For the spin-1/2 Heisenberg Hamiltonian, only six of the sixteen elements of the transfer 
matrix are non-zero, namely those that leave the 3-component of spin unchanged. 
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The 2D Heisenberg model is of interest because it describes the undoped precursor insula- 
tors of high temperature superconductors, and it has been studied with various Monte Carlo 
techniques ||, including a very efficient loop cluster algorithm [H||. The cluster algorithm 
has led to a very precise determination of the low energy parameters of the Heisenberg model 
H, consistent with the chiral perturbation theory treatment of Hasenfratz and Niedermayer 
||. The loop cluster algorithm has a well-defined continuous-time limit, as discussed be- 
low. The computational advantage of operating directly in the time continuum allows us 
to verify the predictions of chiral perturbation theory at temperatures that are inaccessible 
with the conventional discrete-time algorithm. The results for the low energy parameters of 
the theory are consistent with previous calculations as well as with experimental data on 
precursor insulators of high temperature superconductors |7]||. 

The loop cluster algorithm in discrete time is implemented on a space-time lattice with 
(L/a) 2 spatial points and AN time slices, with a ±1/2 spin state defined at each lattice point. 
Both in space and in Euclidean time periodic boundary conditions are used. Decomposing 
the Hamiltonian with the Trotter formula allows us to define an action for each lattice 
configuration. Loops are built up according to a set of rules dictating the loop "flow" through 
each plaquette. After a loop is built, every spin on the loop is flipped. The loop-building 
rules are calibrated to satisfy the condition of detailed balance, ensuring that the path 
configuration space is sampled in proportion to the statistical weight of each configuration. 

The discrete-time rules, together with their continuum analogs, are summarized in Table 
|. Time is in the vertical direction for the plaquettes displayed here. The essence of the 
loop cluster algorithm is that we assign a "flow" to each legal plaquette, that is, to each 
plaquette with a non-zero entry in the transfer matrix and hence with finite action. In the 
case of the "optional decay" plaquettes the flow can proceed in either the space or the time 
direction; we choose the probability for proceeding in the space direction so that the entire 
algorithm satisfies detailed balance. Starting at a random site in space-time, we connect 
flows until the loop closes, and then all spins on the loop are flipped. The cluster-building 
rules ensure that the cluster paths always close. 
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The continuum limit is determined by first observing that jumps from one spin site 
to a neighbor (the "optional" plaquettes in the table) are sporadic and of measure zero, 
so that most plaquettes are continuations in time. The solid lines in Table | represent 
continuously spin up states, and the dotted lines are continuously spin down states. The 
forced transition plaquettes become "bonds" in the continuum limit, and the flow is required 
to travel along the bond and reverse direction in the continuum paradigm. The optional 
decay plaquettes become Poisson-like decays between unlike spins; the probability of a jump 
in the discrete case becomes a probability per unit time in the continuum case. The duration 
of a path segment on a given spin site is thus distributed like the lifetime of a radioactive 
nucleus: the more decay channels available, the more likely the path will terminate in a 
given amount of time. The essential difference between the discrete-time and continuous- 
time implementation is that the former requires us to store spin state information for each of 
the points of the space-time lattice, while the latter requires us only to store the transition 
times for each spin site (plus an extra bit to record the state at t — 0). In the code developed 
here, all the transition times were maintained in a large array. A cyclic double-linked-list 
storage technique was implemented to ensure that the overhead associated with inserting 
and deleting time values was minimized. 

A significant advantage of the loop cluster algorithm is that it allows for the implemen- 
tation of improved estimators. We found that the improved estimators for susceptibility, 
staggered susceptibility, and internal energy density all have easily determined continuum 
limits. 

The continuous-time algorithm with improved estimators was implemented in a FOR- 
TRAN code ||. Instead of the node-by-node method of the discrete-time algorithm, the 
continuous-time algorithm is based on the exponential distribution of decay times. Starting 
at a randomly selected space-time location in the lattice, the neighborhood of the current 
endpoint of the path is surveyed. The first job of the survey is to determine the "hori- 
zon", which is the time of the next forced transition. The other purpose of the survey is 
to determine the number and duration of decay channels — i.e. sites with opposite spin 
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unvisited by the current path. The time of the nearest change in the neighborhood is called 
the "waypoint." A pseudo-random number is generated, and a decay-point is computed ac- 
cording to the exponential distribution corresponding to the number of decay channels. The 
current time marker is moved to the closer of the waypoint or the decay-point. If the marker 
reaches the waypoint, then the waypoint and decay-point determination is repeated. If the 
marker reaches either the horizon or a decay-point, then the path moves to the appropriate 
neighbor, selected randomly from several channels if necessary. The process continues until 
the path closes. 

In general terms, the computer time required to run a simulation for a given volume and 
temperature is roughly comparable for the discrete code (at a usable time-granularity) and 
the continuum code, but the continuum code requires substantially less storage space. The 
real advantage of the continuum code is that it obviates the need to conduct the several 
runs of successively finer time-granularity needed for extrapolation to the continuum limit, 
thus obviating a costly dimension in the simulation procedure. 

In the ground state of the 2D Heisenberg antiferromagnet the staggered magnetization 
develops an expectation value, and hence the 0(3) spin rotational symmetry gets sponta- 
neously broken down to 0(2). The low energy excitations of the systems are spinwaves 
(so-called magnons) which are the Goldstone bosons of the spontaneously broken 0(3) sym- 
metry. Using symmetry arguments, chiral perturbation theory makes very strong predictions 
for magnon dynamics at low temperatures, with three low energy parameters as unknown 
constants. These parameters are the staggered magnetization Ai s , the spin wave velocity 
he, and the spin stiffness p s . The predictions of chiral perturbation theory at extremely low 
temperatures were checked using a range of square volumes with side length L/a — 6, 8 ... 20 
and a range of inverse temperatures (3 J = 1, 2, 5, 10, 20, 30, 40, 50, 80, 100. Note that the very 
small temperatures T ps 0.01J are inaccessible to the discrete-time code, largely due to stor- 
age limitations. In every case 100,000 configurations were generated after a thermalization 
period of 1000 configurations. 

Several of the key predictions of chiral perturbation theory were directly verified with 
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this code. In particular, the energy spectrum is that of an 0(3) rotor, with energy levels 
characterized by an integer spin j, having degeneracy 2j + 1 and energy levels distributed 
as j(j + 1). Finite volume effects are computed as expansions in he/ p s L, in units where the 
lattice spacing a = 1 . The energy spectrum is computed in reference M to be 



Ej — E = 

(hA 2 T hr r Qnn9fiK 

(7) 



2p s L 2 



he 3.900265 „, r 

+ 0(L- 2 ) 



p s L Aix 

Including only the leading term in the partition function, the uniform susceptibility ap- 
proaches the form 

and the staggered susceptibility goes to the temperature-independent form 



t^o (he) 2 



v , ^HlH^l T 4 

Xs =— t JZ xo ^ 



3.900265 „, r 2 , 
1 + 3— + 0{L~ 2 ) 



(9) 



These functional forms in the large (3 limit, including both volume and temperature de- 
pendence, were verified for uniform and staggered susceptibility. In addition, the staggered 
susceptibility is shown to plateau at large (3 as predicted by chiral perturbation theory. 

An independent fit for the chiral perturbation theory parameters M. S) he, and p s gives 
excellent agreement with the results of 0, as shown in Table ||. This fit required that the 
partition function of the 0(3) rotor be included in its entirety, instead of including just the 
leading term, as the limiting forms in equations (§) and (0) employ. In addition, only inverse 
temperatures (3 J > 10 were used in the fit, since the rotor approximation is valid only for 
very small temperatures. In order to reproduce the accuracy of the fit in reference 0, it was 
necessary to find fitted values for the coefficients of the quadratic terms in the expressions 
for energy and staggered susceptibility, that is, the coefficients of the terms (hc/p s L) 2 in 
the expressions above. Note that the agreement between the current work and reference 
is particularly remarkable because they are derived for different volume-temperature 
regimes. Reference || was concerned with the "cubical" regime TL/hc = 1, while the 



current study focuses on the "cylindrical" regime TL/Tic « 1. Quadratic coefficients for 
the cubical regime are not included in Table [TT] because they are not comparable to those 
in the cylindrical regime. The five parameter values in Table [TT] resulted in a goodness-of-fit 
X 2 /d.o.f. = 1.5. Figures |l](a) and 0(b) show the fit for uniform and staggered susceptibility, 
respectively. Solid lines represent the fitted chiral perturbation theory result. Circles and 
error bars are displayed at each simulation point. The fit is quite good for (3 J > 10. 

In conclusion, we find that the continuous-time formulation discussed here provides a 
superior method for evaluating path integrals that arise in the study of discrete quantum 
systems. Using this method allows us to validate the chiral perturbation theory results 
for the Heisenberg antiferromagnet in the extreme low temperature regime. In discrete 
Euclidean time a loop cluster algorithm has also been constructed for lattice fermions ||10|| . 
Our continuous time method can be applied to this problem as well. In fact, in principle 
it works for any discrete quantum system, even if a cluster algorithm may not be available. 
For example, for a nonrelativistic particle living in a finite volume one can formulate the 
path integral using the discrete momentum states. Similarly, for a lattice gauge theory with 
a compact gauge group one can work in the discrete basis of representations. In both cases 
this yields a well-defined path integral directly in the Euclidean time continuum. 
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Alfred P. Sloan Fellowship (UJW). This work has also been supported in part by funds 
provided by the U.S. Department of Energy (D.O.E.) under cooperative research agreement 
DE-FC02-94ER40818. 
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FIGURES 

FIG. 1. (a) Uniform susceptibility and (b) staggered susceptibility versus inverse temperature 
for various volumes. Solid lines are predictions from chiral perturbation theory with fitted param- 
eters; circles with error bars are lattice simulations with continuous time. The fit is very good for 
(3 J > 10. 
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TABLES 

TABLE I. Summary of Plaquette Flow Rules. +1/2 and solid lines indicate spin-up sites; —1/2 
and dotted lines are spin-down. The time direction is vertical. Flow rules for inverse plaquettes 
are analogous. 
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TABLE II. Comparison of Fitted Parameters 



Ref § Current 

Spin stiffness p s 0.186(4) 0.185(2) 

Spin wave velocity he 1.68(1) 1.68(1) 

Staggered Magnetization M s 0.3074(4) 0.3083(2) 

Quadratic coeff.: Energy 0.068(1) 

Quadratic coeff.: Stag. Susc. 0.338(7) 
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This is Figure 1 (a) for the preprint CTP-2509, "Simulations of Discrete Quantum Systems in Continuous Euclidean Time" by B.B. Beard and U.-J. Wiese 
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FIG. 1 (a) Uniform susceptibility versus inverse temperature 
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perturbation theory with fitted parameters; circles with error 
bars are lattice simulations with continuous time. The fit is 
very good for (3 J >10. 
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This is Figure 1 (b) for the preprint CTP-2509, "Simulations of Discrete Quantum Systems in Continuous Euclidean Time" by B.B. Beard and U.-J. Wiese 
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FIG. 1(b) Staggered susceptibility versus inverse 
temperature for various volumes. Solid lines are 
predictions from chiral perturbation theory with 
fitted parameters; circles with error bars are lattice 
simulations with continuous time. The fit is very 
good for (3J >10. 
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